Modulation instability and solitary wave formation in two-component Bose-Einstein 

condensates 
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We investigate nonlinear dynamics induced by the modulation instability of a two-component 
mixture in an atomic Bose-Einstein condensate. The nonlinear dynamics is examined using numeri- 
cal simulations of the time-dependent coupled Gross-Pitaevskii equations. The unstable modulation 
grows from initially miscible condensates into various types of vector solitary waves, depending on 
the combinations of the sign of the coupling constants (intracomponent and intercomponent). We 
discuss the detailed features of the modulation instability, dynamics of solitary wave formation, 
and an analogy with the collapsing dynamics in a single-component condensate with attractive 
interactions. 

PACS numbers: 03.75.Mn, 03.75.Lm 



I. INTRODUCTION 

Dynamics of spatial pattern formation in nonlinear me- 
dia is important for a wide range of physical phenom- 
ena 0. Modulation instability (MI) is an indispensable 
mechanism for understanding pattern formation from a 
uniform medium. MI occurs when a constant-wave back- 
ground becomes unstable to induced sinusoidal modu- 
lations under the combined effects of nonlinearity and 
diffraction in a spatially nonlinear field. Also, it is known 
that MI causes a uniform medium to break-up into pulsed 
"solitary waves" 0. The effect of self-interactions in 
nonlinear media plays a crucial role for the MI in a sin- 
gle component system. In a multicomponent system, in 
which there is more than one order parameter, additional 
types of interactions can occur between different compo- 
nents and have a great influence on the MI. 

Atomic-gas Bose-Einstein condensates (BECs) are a 
good system for examining MI and nonlinear matter- 
wave dynamics. In this system, the nonlinearity origi- 
nates from the atom-atom interaction. Moreover, manip- 
ulation of the matter waves can be achieved by applying 
established techniques from atomic, molecular, and opti- 
cal physics. For example, in single-component BECs, re- 
searchers have studied nonlinear excitations such as dark 
solitons || p, P , bright solitons j(| 0, Q , and quantized 
vortices jjj [Hj • It has been clarified that the MI plays 
a crucial role in the formation process of those excita- 
tions a nd g ives rise to the intriguing nonlinear dynamics 

This paper addresses the nonlinear dynamics of soli- 
tary wave formation induced by MI in a mixture of two- 
component atomic-gas BECs. The experimental creation 
of multicomponent BECs has been achieved by the si- 
multaneous magnetic trapping of the atoms lying in the 
week- field- seeking states 0, ^J, the use of an opti- 
cal trapping that liberates the spin degrees of freedom 
of atoms | l8l 1 1 9| or simultaneous trapping of different 
species of atoms [2(3, El ■ The MI in two-component 
BECs was firstly discussed by Goldstein and Meystre 



[23| . and recently reexamined |24|. However, how the 
condensates develop under the MI is still unclear, though 
some progress can be seen for the study of spin-1 BECs 

H 

First, we summarize the MI condition with respect 
to two intracomponent interactions and an intercompo- 
nent one. Although all of these atom-atom interactions 
were repulsive in the past experiments of multicompo- 
nent BECs, except for a boson-fermion mixture in Ref. 
|22| . further insight can be gained if we consider some 
of them are attractive. The character of the interac- 
tions ma y b e controlled by choosing the specific kinds 
of atom pel l2lL |22| or by using a homonuclear ||| Q 
or heteronuclear Fcshbach resonance prl |27|. Next, we 
discuss the nonlinear dynamics caused by the MI, empha- 
sizing the role of the intercomponent interaction. We fo- 
cus on the situation in which one component is suddenly 
put on the other component with repulsive interaction. 
The MI will lead to the formation of a vector soliton 
train in such a way that a bright soliton train is gener- 
ated thr oug h the MI in a single-component condensate 
[Til ll2L Il3j . Our previous paper [2g revealed the dy- 
namics of domain formation of two-component BECs in 
the case of all interactions being repulsive, in excellently 
agreement with the experimental observation of Miesner 
et al. [23 . We extend the analysis to the cases where two 
components have attractive intercomponent interaction, 
or one of them has an attractive intracomponent inter- 
action. Depending on the combination of the sign of the 
s-wave scattering lengths, two-component BECs exhibit 
rich nonlinear dynamics of solitary wave formation. 

This paper is organized as follows. In Sec. [n] we for- 
mulate of the problem for two-component BECs using 
a quasi-one-dimensional model for simplicity. Then we 
use linear stability analysis to clarify the MI with re- 
spect to the sign of the coupling constants. Section ITTT1 
presents the numerical simulation results that confirm 
the MI analysis and show how the solitary wave forma- 
tion occurs in the condensates through the MI. Section 
IIVI is devoted to conclusion. 
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II. MODULATION INSTABILITY OF 
TWO-COMPONENT BOSE EINSTEIN 
CONDENSATES 

A. Model 

We start with a two-component BEC with atomic 
masses mi and m 2 . The dynamics can be derived by 
assuming that the two condensates are described by the 
wave functions ^i(r, t) and ^2 (r, t) . At zero tempera- 
ture, the total energy functional of the system is 




+^|*i| 2 + KSl^l 2 + ^il^il 4 + ^l^l 4 

+ <?12|*l| 2 |* 2 | 2 • (1) 

The condensates are assumed to be trapped in axisym- 
metric harmonic potentials: 

V c ^(r,z)^^ mi uf(r 2 +\ 2 z 2 ), » = 1,2, (2) 

where u>i is the transverse trapping frequency and A is the 
aspect ratio of the potential. Each component can have 
its own trapping frequency due to the ^-factor and index 
of the atomic hyperfine levels along the quantized axis. 
The intracomponent coupling constant gi = 4Trh 2 ai/mi 
is characterized by the scattering lengths a\ and a 2 be- 
tween atoms of the same species, while the intercompo- 
nent one g\2 = 'inh 2 0,12 / 'm.12 (m^ 1 — rn^ 1 + m^ 1 ) is 
determined by the scattering length a\2 where an atom 
in the ^1 component scatters from another atom in the 
^>2 component. This intercomponent coupling yields new 
structures and dynamics not found in a sin gle c omponent 

bec EE E3, HE IM m, in, [33, H3, [M lali . 

The dynamics of two-component BECs can be de- 
scribed using the coupled GP equations, which are de- 
rived from the variational principle ihd^i/dt = 8E/5^* 
as 

lh ^f = (-^ + ^+»l*a| a +5i2|*i| a )*d(3b) 

The normalization of each wave function is taken inde- 
pendently as J dr|\I>i(r)| 2 = iV,-. 

We assume that the condensates are tightly confined in 
the transverse direction, so A <C 1. This condition means 
that the motional degrees of freedom in the x-y plane 
are frozen, a situation that could be realized in highly 
elongated cigar-shaped potentials. In this case, one can 
factorize the condensate wave function into a longitudinal 
and a transverse part as 

* i (r,t)=4> ( i\x,y)ip i (z,t)e- i "* t , (4) 



where <p\_ (x, y) is the normalized ground state of the 

transverse potential v}^(r) = rriiLufr 2 /2 with energy 
hu>i. The system is thus effectively reduced to a one- 
dimensional geometry, with the longitudinal condensate 
wave function ipi{z,t) satisfying the one-dimensional GP 
equations: 

a^i / h 2 d 2 1 , 2 2 2 , I , ,2 

+ui2|^2| 2 ^V>i, (5a) 
dip 2 ( h 2 d 2 1 222, 1 , ,2 

+"i2^i| 2 ^2- (5b) 

Here, 

u t = g t m = g t J dxdy\4> ( ^\ A = (6) 
U12 = .9i2?/i2 = .912 J dxdylcj)^] 2 ]^] 2 = ^fq: jpj ( 7 ) 
with the length scale hi = ^h/rriiUJi characteristic of 

B. Relation with the nonlinear Schrodinger 
equation in nonlinear optics 

In the context of nonlinear optics, a special attention 
has been paid to MI in Kerr media in which light-wave 
propagation is described by the nonlinear Schrodinger 
equation (NLSE) within the scalar approximation of the 
electromagnetic field 2]. The NLSE exhibits instability 
of self-phase-modulation (SPM) when nonlinearity and 
group velocity dispersion (GVD) act in opposition, e.g., 
for self-focusing waves associated with negative nonlin- 
earity the GVD should be "normal" (a positive GVD 
coefficient) and for self-defocusing waves associated with 
positive nonlinearity the GVD should be "anomalous" (a 
negative GVD coefficient). This condition is also neces- 
sary for the existence of bright solitons which result from 
an exact balance between nonlinearity and dispersion. 

If accounting for polarization of the electromagnetic 
field, light propagation in isotropic Kerr media is de- 
scribed by two incoherently coupled NLSEs instead of 
the single NLSE |2J. Then, the incoherent coupling be- 
tween two NLSEs, referred to as cross-phase-modulation 
(XPM), leads to MI for any sign of nonlinearity and GVD 
|34| . The XPM is a general phenomenon characteristic of 
the simultaneous nonlinear propagation of several waves 
belonging to different modes. Also, MI induced by the 
XPM is of fundamental importance as it suggests the 
possibility of soliton formation in the normal dispersion 
regime. 
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Equations (0 have a close analogy with the incoher- 
ently coupled NLSEs in nonlinear optics, where the itj- 
and Mi2-terms correspond to the SPM and XPM terms, 
respectively. In nonlinear optics, the ratio of the non- 
linear coefficients for SPM and XPM can be altered us- 
ing the light's angle of elliptic polarization Q. For the 
atomic BECs, the strength of the atomic interactions can 
be altered using the Feshbach resonance [26|, |27j ■ 



C. Modulation instability analysis 

In a single-component nonlinear wave, the MI induced 
by SPM exists only for the waves with self-focusing non- 
linearity, corresponding to the attractive interaction be- 
tween atoms. The intercomponent coupling (i.e., XPA) 
is a feature of the two-component system that does not 
exist in a single-component system. In this section, we 
discuss how the sign and strength of the coupling param- 
eters u\, u 2 and iti 2 affect the MI. 

We examine the stability of miscible two-component 
BECs with the homogeneous one-dimensional density 
mo = |-0io | 2 and n 20 = 1-020 1 2 0, IF 
wave functions are written as tpi(z, t) - 
the linearized equation for the fluctuations becomes 



2g. When the 

/r^ + 5ipi(z,t), 



at 



at 



h 2 d 2 



Sip! + umioWi + Sipl) 



2m i dz 2 

+«12 V n 10^20 {Sip2 + H2) ( 8a ) 



h 2 d 2 



5ip 2 + U2n 2 o(Slp2 + S1P2) 



r\ 7 9 " r ^ 1 "/'"ZUV 

, ~^2) 

2rn 2 dz z 

+ui2^n w n 2 o(Sipi + Sipl) (8b) 



We assume a general solution of the form Sipi — 
(i cos(kiZ — Sit) + irji sin(fciZ — Sit), where we allow for 
different wave numbers fcj for the ipi (i = 1, 2) compo- 
nents. Then, Eqs. JSJ provide a set of equations for the 
amplitude Q and Straightforward calculation gives 
the dispersion relation 



(SI 2 - A^iSl 2 - A 2 ) = P 2 



where 



P = 



2 1,2 



h 2 k, 

2m 5 ; ^ 2m 
U12 



2inn 



V"l0«20^lfc2- 



(9) 

(10) 
(11) 



The dispersion relation gives a quadratic algebraic equa- 
tion in terms of SI 2 , whose solution is 



St 



1 



Ax + A 2 ± v /(Ai + A 2 ) 2 + 4(P 2 -A 1 A 2 ) 



The condensates are uniformly miscible and their sta- 
bility is governed by Eq. (|12|) . If the frequency Sl± has 
an imaginary part, the spatially modulated perturbations 
grow exponentially with time, as is evident from the form 



of Sifti. This unstable growth of weak perturbations is re- 
ferred to as the MI. The MI condition depends on the 
sign of two variables 

A = Ai + A 2 , (13) 
A = P 2 -AiA 2 ; (14) 

With these two variables, Eq. i|12fl is rewritten as 

Sl\ = ~(A± 0Y 2 + 4A). (15) 

For A > 0, the value of Sl 2 + is always positive, whereas 
the Sl 2 _ becomes negative only if A > in which case Sl_ 
is purely imaginary. For A < 0, the value of Sl 2 ^ becomes 
negative when A < and Sl 2 _ is always negative; thus, 
the system is always modulationally unstable. 



D. Condition of MI for a single-component 
condensate 



Before considering the general case, it is instructive to 
review briefly the MI of a single-component BEC. For 
a single component, u 2 = 0, ui 2 = (P = 0), and 
the dispersion relation Eq. I|12|l reduces to SI 2 = Aj. 
The MI occurs when Ai = h 2 k 2 /2m\ + 2uin w < 0. 
We obtain the well known result that the MI occurs 
only with an attractive coupling constant u\ < 0. In 
this case, the imaginary component of the frequency 
G = ImSl represents the growth rate of the modula- 
tion, which is called the gain spectra This compo- 
nent is given by G — -^^^ ^ m i\ u i\ n w — Ti 2 k 2 in the 

range < k\ < ,/4m~iJuiJnw/h. The fastest growth 
occurs for the wave number fci max that gives a maxi- 
mum of G. The extremum condition dSl 2 /dk 2 = gives 
fcimax = \/2mi\ui\nio/h and the maximum growth rate 
G m ax = \ui\nio/h. The MI associated with the attrac- 
tive interaction has a key role in the formation of bright 
solitons of a sing le-component BEC [UHl. 



E. Condition of MI for a two-component 
condensate 



This study is concerned with the MI relevant to the 
intercomponent coupling; thus, we fix the interaction 
of the ^-component to be positive u x > 0. Possible 
choices of the sign of the coupling strengths w 2 and u\2 
are summarized in Table [I] To classify the types of the 
instability more clearly, we introduce the length scale 
i 2 = (4toiui?iio / fi 2 ) -1 and the dimensionless wave num- 
ber fcj = k t i. Then, Eqs. ljT3|) and ((HI) become 



A = 



~k 2 (k 2 + l) + ^~k 2 (k 2 + l2 ) 



A = 



/ h 2 



V4m 2 ^ 4 



k 1 k 2 



7i2 



(k 2 + l)(k 2 



■72) 



(16) 
(17) 
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TABLE I: Four cases that are considered, each defined by the 
sign (positive = "+h, negative = "— h) of the coupling con- 
stants U2 and iti2. In all cases, ui is assumed to be positive. 

112(72) ^12(712) 



(a) 
(b) 
(c) 
(d) 



+ 
+ 



+ 



+ 




FIG. 1: The modulationally unstable region in fci-fe space 
for the combinations (a) and (b) in Table Q The un- 
stable region lies below the thick line boundary given by 

h = V / 7i 2 2 /(^ + l)-72- The magnitude of the gain spectra 
G = Imf2 (arbitrary unit) is shown by a contour plot in the 
unstable region. 



where 



72 



m 2 u 2 n2o 
miUiriio 



712 



"12 
Ml 



m 2 n 2 o 
minio ' 



(18) 



These equations show that the MI condition depends on 
two atomic masses, condensate density, three coupling 
constants, and the range of the wave numbers. Here, we 
assume mi = m 2 , which greatly simplifies the form of 
the following equations. 

We search the unstable region of the wave number fc; 
by changing the values of 72 and 712. By examing the 
possible choices for the signs of 72 and 712 shown in Table 
[I] we obtained the unstable region in ki~k 2 space as shown 
in Fig. 2] and [21 We summarize below some features of 
the instability. 



Region (a): 72 > and 712 > 



This situation is of particular importance because the 
individual condensates with the repulsive interaction are 
modulationally stable. In this case, A > and thus 
f2j_ is always positive. Therefore, the unstable condi- 
tion is determined from Eq. (|17f) . where the positive 
A (A > 0) gives a purely imaginary f2_. For fixed 
k 2 the instability occurs within a certain range of k\; 





< fci < yli 2 /(k 2 +72) — 1 as shown in Fig ^ F° r 
the wavenumber to be real, the term in the square root 
must be positive, which gives the necessary condition 

72 



712 > \/72 + k 2 01712 < -^72 + k\ for the MI to occur. 
The former corresponds to the strong repulsive intercom- 
ponent interaction and the latter to the corresponding 
attractive interactions. Thus, the unstable range is inde- 
pendent of the sign of 712. In the former case, we obtain 
the well-known condition ^/gig 2 < g\ 2 for phase separa- 
tion in the long wavelength limit ki — > [3(|. In Fig. ^ 
we also show the magnitude of the imaginaiy component 
of f2_ (gain spectia G = Imf2_). The maximum of G ap- 
pears at the wave number k\ =k 2 = fc max . After setting 
k\ = k 2 = k in Eq. I|15|l . the most unstable wave numbei 
is calculated from dVl 2 _/dk — 0, with the result 



(72 - I) 2 + 4 7l 2 2 - 72 -1 



1/2 



(19) 



and the maximum growth rate becomes 



hk 2 



2ml 2 8m£ 2 



V (72-l) 2 



47 



72 



m 

The unstable modulation develops by following the 
eigenvectors associated with the eigenvalue For 
ki = k 2 = k, a simple calculation gives the mode am- 
plitudes as 




^± signal + ^ 



(21) 



Vi± 



1 



v± 



1,2. (22) 



The positive (negative) sign represents the eigenvector 
associated with Q + (^-)- F° r positive u\ 2 the amplitude 
Ci+ (Ci-) i s always positive (negative), which means that 
the unstable modulation develops into out-of-phasc 
waves. This feature follows from the fact that the repul- 
sive character of the intercomponent interaction forces 
the two components apart. 



2. Region (b): 72 > and 712 < 

In this case, although each component has a repulsive 
intracomponent interaction, the two components have a 
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FIG. 2: The modulationally unstable region in fci-fe space for 
cases (c) and (d) in TableU The unstable region lies below the 

thick line given by k<2 = y 712/(^1 + 1) — 72- The boundary 
curve crossing the origin (bottom left) is given by Ai + A2 = 0; 
the right (left) region from the boundary represents the region 
(i) ((ii)) ( see t ne text)- The contour plot shows the magnitude 
of the gain spectra G = Imf2 in arbitrary unit. 

strong attraction. Because the combination of the inter- 
component coupling is included through 712, this case is 
similar to that of case (a). Thus, the most unstable wave 
number and the corresponding maximum gain spectra 
are the same as the combination (a) with only the signs 
of the modulation amplitudes being different: Q+ < 
and £,-_ > in Eq. J2IJ- Therefore, the MI leads to an 
in-phase evolution of the two-component modulation. 

3. Region (c): 72 < and 712 > 

Next, we consider the situation in which one compo- 
nent has an attractive intracomponent interaction. As 
shown in Sec. Ill Dl the ^-component always undergoes 
the MI induced by the attractive force. It is interest- 
ing to consider how the presence of the other component 
affects the MI condition. The unstable region in k\-k2 
space and the gain spectra G = ImSl is shown in Fig. [21 
For 72 < there appears to be a region satisfying A < 0, 
where f2+ can become negative (i.e., unstable) if A < 0. 
However, there is no region where these two inequalities 
(A < 0, A < 0) are satisfied simultaneously. Hence, we 
will focus on the instability associated with Q 2 _ . 

For 712 = 0, independent of the values of fc 1; there is 
an "unstable band" with G = ^72 — ^2 m the range 
< ft2 < VTY^T- For 712 7^ the unstable region extends 
beyond the boundary line ft2 — -v/ 1 T2 ] - Two unstable 
regions exist: (i) A > with A > and (ii) A < 



with A > 0. The condition for region (i) is the same 
condition as that for case (a) except the boundary line 
k,2 = \J1\2I (^i + 1) + I72I has no intersection with the 
fci-axis because 72 < 0. The condition for (ii) yields the 
imaginary f2_ in the region bound by the curves k\ = 
and Ai + A 2 = in Fig. [3 which is comprised by the 
region given by (i). As a result, the MI condition is solely 
determined by A > as in case (a). 

In region (c) , the system undergoes MI induced by both 
the intracomponent interaction U2 and the intercompo- 
nent one U\2- The most unstable wave number and the 
maximum growth rate are again given by Eqs. (|19|l and 
(J2DJ. When 712 = 0, we reproduce the results of the 
single-component case G max = ?i|72|/4mi? = \u2\ri2o/h. 
Therefore, an increase of 712 always increases the un- 
stable growth rate over that of the single-component 
case. When G m ax is located within the unstable band 
< h,2 < a/|72 I, the dominant contribution to the MI 
should be the intracomponent attraction associated with 
the negative U2- If 712 > I72K2I72I + 1), the location of 
G max moves outside the unstable band, so that the inter- 
component repulsion it 12 has the dominant influence on 
the MI. Compared with the case (a), the magnitude of 
the gain spectra is larger by a few factor of about 0(1). 
This is due to the multiplication effect of those two in- 
stabilities. 

Once the instability starts, the modulation develops by 
following the eigenvectors associated with the eigenvalue 
fi_. According to Eq. I|21|l . the modulation becomes 
out-of-phase. 

4- Region (d): 72 < and 712 < 

This region is similar to that of (c). The MI condition, 
the most unstable wave number and the corresponding 
maximum gain spectra are the same as those in region 
(c). From the eigenvectors of the modulation amplitudes, 
the MI is related to an in-phase evolution of the two- 
component modulation. 



III. NUMERICAL SIMULATIONS 
A. Formulation of the simulations 

In this section, we present the results of our numerical 
simulations that illustrate the effect of the MI on the 
nonlinear evolution of the condensates. First, we describe 
the formulation, the initial conditions and the parameters 
of the simulations in detail. 



1. The dimensionless GP equations with particle loss terms 

To reduce the number of the parameters, we assume 
mi = m% = m and lo\ — L02 = so that 61 = &2 = 



6 



&ho = \/h/mLUj_. It is also convenient to introduce the 
scales characterizing the trapping potential wj 1 , &i 10 and 
hu>± for time, length and energy, respectively. By replac- 
ing the wave function with the total particle number N 
(= Nx + N2) as ipi —> ipiy/N/bho, the one-dimensional 
GP equations © reduce to 



Ml 


■ 1 d 2 




~2d^ 


.di>i 


" 1 d 2 




~2~d^ 



X 2 

+ ^z 2 + u 1 \4> 1 \ 2 + u 12 \i)2\ 2 

X 2 

+ yZ 2 + C/ 2 |V 2 | 2 + C/l2|V'i| S 



V-i(23a) 
^#3b) 



where Ui — UiN /hio±b\ w = 2Nai/b\ 10 (i — 1,2), C/ 12 = 
U12N ' 1 'fkv j_&ho = 2iVoi2/6h 0) and the normalization of the 
wave function is J dz\tpi(z)\ 2 = Ni/N. 

We will simulate the dynamics for the condensate with 
an attractive interaction. Therefore, we should include 
an effect of the atomic loss due to inelastic collisions [^J ■ 
We model this effect by adding to Eqs. (|2*3*)l the phe- 
nomenological loss term 



loss term = —i ( Lj \tpj\ 4 



3-jl 



$3 



with 



1 T^N 2 

f(3) = 1 L 3 iV 

s 3!6tt 2 u;±&£ ' 



Ldn — 



1 L dii N 
2! 4ttc^6? 10 



J = 1,2 
(24) 

(25) 



in the right-hand side of Eq. (|23[) . The first term on 
the right side of Eq. (|24|l is related with the three-body 
inelastic collisions, which is the dominant mechanism of 
particle loss when the self-focusing collapse of the con- 
densate occurs. The second term on the right side of 
Eq. I|24|) represents the inelastic loss due to collisions 
between different components, which is associated with 
inelastic collision between different atomic species |3(| or 
the spin exchange collision for the two components with 
different hyperfine states [33 • Because the detail of the 
particle loss through the collapse are not needed here, 
we use a value for Lj and Ldif that is consistent with 
experimental results: l£_ = 1 x 10~ 26 cm 6 /s [38| and 
L dif = 3 x 1CT 14 cm 3 /s j^. 



2. The initial conditions 

We numerically solved the time-dependent GP equa- 
tions H23fl using a Crank-Nicolson implicit scheme with 
8 x 10 3 grid points and a time grid At = 5.0 x 1CT 4 . The 
focus here is on how the MI grows spontaneously from the 
miscible condensates. Therefore, the initial two compo- 
nents should be uniform and overlap as much as possible 
in the trapping potential. To do this, we first prepared 
the stationary solution of component (denoted by ^i n i) 
by propagating Eq. (|23b 'l in imaginary time, under the 
normalization J dz\tpi\ = 1. This was done without the 
[/12-term and the particle loss term. Next, at t = 



in real time simulations, some fractions of ip\ compo- 
nent were suddenly put into the f/ ; 2 component that had 
the same density profile as tf>i. Initially, each compo- 
nent has the same inverted-parabola density profile, but 
a different normalization condition J dz\ipi\ 2 = Ni/N . If 
Ui = U2 = U12 is not satisfied, this initial configuration 
is nonstationary and thus can develop following the MI. 
Hence, we focus on such nonstationary cases. This situa- 
tion can be realized experimentally by using a rf-pulse to 
transfer the population from one hyperfine level of atoms 
to the other one [IE El. 



Parameters and the validity of the one-dimensional 
simulations 



We consider the quasi-one-dimensional geometry char- 
acterized by the aspect ratio A = 0.02. Using the mass 
of rubidium atoms and the radial trapping frequency 
■jj±_ = 2ir x 100 Hz, we obtain the length scale b^ = 1.1 
/im and the time scale luJ 1 = 1.59 msec. According to 
the values of typical alkali atoms, we fix the intracompo- 
nent s-wave scattering lengths for the simulations of the 
cases (a) and (b) in Table [Q as 



<Zi = 5.5 nm, 02 = 5.8 nm 
and for the cases (c) and (d) as 

5.5 nm, 02 = —0.2 nm. 



(26) 



(27) 



Further simplification can be obtained if we confine our- 
selves to distribute the equal particle number for two 
components Ni = N 2 , i.e., J dz\ipi\ 2 — 1/2 (i — 1,2). 
Thus, we have the total particle number N and the value 
of the intercomponent s-wave scattering length 012 as 
variable parameters. 

We used the quasi-one-dimensional model under the 
assumption that the transeverse motion would be frozen. 
To justify this assumption, the energy scale of the trans- 
verse confinement should be much larger than the non- 
linear interaction energy. This yields the condition 
diN/bho *C 1 Hj|. Unfortunately, this condition is not 
satisfied for our parameters; we obtain ciiN /b\ w ~ 1 for 
the typical parameters presented below. However, the 
resulting transverse motion is only a rapid breathing os- 
cillation, which does not affect the Mi-induced dynam- 
ics in the longitudinal direction (28|. The more critical 
condition for our study is to prevent the transverse col- 
lapse 0, being given by | SiraiN \ipi{z)\ 2 /&ho| < H.7. In 
our situation, this position-dependent condition is nearly 
satisfied whenever the focusing collapse of the attractive 
condensates occurs at the trapping center in the following 
discussion. 



B. Numerical results for region (a) 

We first consider region (a) in which all three scat- 
tering lengths are positive. The initial state is given by 
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FIG. 3: (Color online.) Numerical results of Region (a) for 
N = 5 x 10 3 . (a) Time development of condensate density 
|?/>i| 2 (upper panel) and \ip2\ 2 (lower panel) at z — for ai2 = 
5.6 nm, 5.8 nm, and 6.0 nm. (b) The contour plots of the 
density of both components with respect to time ([0,1200], 
horizontal axis) and z ([-100,100], vertical axis) for a\2 = 
6.0 nm. (c) The density profiles of |^i| 2 (solid-curve), IV^I 2 
(dashed-curve) and total density tit = \tpi\ 2 + \ip2\ 2 (dotted- 
curve) for 012 = 6.0 nm at t = 240, t = 580, and t = 960. 
The unit of time is UiT 1 . 



the miscible condensates with the same density profile of 
inverted parabolas. This region, together with the ini- 
tial state, coincides exactly with the experiment of Mies- 
ner et al. |2fl| . They first prepared all 23 Na atoms in 
the \F = 1, to.f = 1) hyperfine state in an optical trap 
and then placed instantaneously half of them into the 
\F = l,mp = 0) state using an rf field. Letting the sys- 
tem evolve freely while using the quadratic Zeeman effect 
to prevent the \F = l,rriF = —1) component from ap- 
pearing, they found that spin domains formed with two 
components alternatively aligned from the initially mis- 
cible condensates. Our previous paper |28| pointed out 
that this observation is due to the MI caused by the in- 
tercomponent repulsive interaction. In the following, we 
describe the features of the nonlinear dynamics in more 
detail using one-dimensional simulations. 



1. Dynamical features 

The MI changes greatly the behavior of nonstation- 
ary development of the condensates. In Fig. |3 we show 



the results of the numerical simulation for N = 5 x 10 3 . 
Figure Efa) represents the development of the conden- 
sate density at the center (z — 0) for several values of 
a\2- A crucial difference of the dynamical behavior is 
seen across the critical value about a\ 2 = \Ja-\a<2 = 5.65 
nm, which corresponds to the criterion for phase sepa- 
ration. When di2 is smaller than the critical value 012, 
|V>i(0, t)\ 2 (1^2(0, t)\ 2 ) first increases (decreases) gradu- 
ally and makes a slow oscillation. Because a 2 > a± and 
O12 > 0, the density of the ipi component is located at 
the center, surrounded by that of the extended ip 2 com- 
ponent. As 012 increases, the oscillation becomes non- 
periodic; the amplitude of \ipi\ 2 drops to zero after some 
time [solid and dashed curves in Fig. a)], which rep- 
resents that ipi is replaced by tp 2 at the center. This 
replacement shows the onset of the ML 

The dynamical process of the spatial pattern forma- 
tion induced by MI is clearly seen in the evolution of the 
overall density. Figure|31(b) and (c) show the evolution of 
each condensate density for 012 =6.0 nm exceeding the 
critical value. Throughout the dynamics, the modulation 
of the density causes two components to become out-of- 
phase. Hence, the total density tit = \ipi\ 2 + \ip2 1 2 keeps 
approximately its initial shape in spite of the irregular 
profile of each component, as shown in Fig. Etc). Af- 
ter t = 300 the density breaks up into smaller domains. 
The domains of the two components alternate in loca- 
tion, while the total density hardly changes even after 
the domain formation. Because a\ < the occupied 
region of \ip2\ 2 expands rather than |'0i| 2 as seen in Fig. 
OUb) and (c). The inelastic loss shrinks the size of both 
components monotonically. 

When the initial particle number increased, the dy- 
namics become more dramatic. Figure Sta) shows the 
time evolution of the central density for N — 5 x 10 4 . 
Compared with Fig. Ola), the central density makes a 
more rapid and complex oscillation after the MI occurs. 
Then, the condensates form much more domains than 
those for N = 5 x 10 3 . Figures Efb) and (c) show that 
the amplitude of the density modulations grows first near 
the edge of the condensate. This growth proceeds from 
the edge to the center, leading to the formation of local- 
ized condensate domains. Since the total number N is 
large, the inelastic loss shrinks the condensate size much 
faster than the case of smaller N. 

Spatially localized domains can be created even in con- 
densates with a repulsive interaction. This is a salient 
feature in a multicomponent system; in the case of a 
scalar condensate (without a periodic potential), local- 
ized domains such as bright solitons are created only 
when the interaction nonlinearity is attractive. Actually, 
the generated domains have a solitary wave structure, 
where the spatial distribution of the condensate phase 
9i =axgipi is almost flat within each domain and its value 
jumps across the density dips where the domains of the 
other exist ■ 

After multiple domains form, the dynamics of each do- 
main may be determined by the phase-dependent inter- 
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FIG. 4: (Color Online.) Numerical results of region (a) for 
N = 5 x 10 4 . (a) Time development of condensate density 
\ipi\ 2 (upper panel) and \ip2\ 2 (lower panel) at z — for N = 
5 x 10 4 and ai2 = 5.6 nm, 5.8 nm, and 6.0 nm. (b) Contour 
plots of the density of both components with respect to time 
([0,1200], horizontal axis) and z ([-160,160], vertical axis) for 
a\2 = 6.0 nm. (c) The density profiles of |i/>i| 2 (solid-curve), 
\1p2f 2 (dashed-curve) and total density nr = \tpi\ 2 + \ip2\ 2 
(dotted-curve) for ai2 = 6.0 nm at t = 200, t = 600, and 
t = 1000. The unit of time is uT 1 . 



action between intracomponent domains and the density- 
dependent interaction between intercomponent domains. 
It is known that the interaction between the bright soli- 
tons in a single component BEC depends on their phase 
difference hTlhJ. In our simulation, the evolution of the 
phase difference between domains is determined nontriv- 
ially, following the nonlinear dynamics caused by the MI. 
When two domains of the same component approach and 
share the same spatial location, the domains exchange 
particles. However, the domain of the other component 
blocks this approach because of the repulsive mean-field 
interaction between domains of different components. 
This is a phenomena analogous to the Josephson effect, 
predicted in Ref. [4fJ, where two single-component con- 
densates with the different phases are separated by a po- 
tential barrier. The oscillations in Figs. and Rafter the 
MI occurred may be caused by the cooperative oscillation 
of the two-component soliton trains by this Josephson ef- 
fect, but this needs further investigation. 



2. Comparison with the MI analysis 



The above dynamics were triggered by the MI induced 
by the intercomponent coupling U\ 2 oc a\ 2 . The first 
modulation growth is determined by the fastest growth 
mode that has the most negative value of O 2 . of Eq. Ijl5|l • 
The corresponding wave number fc max and the maximum 
gain G max will determine the early behavior of the dy- 
namics such as the modulation growth time and the num- 
ber of initially created domains. If we assume the homo- 
geneous condensates, they are given by Eqs. l(T§j l and 
()20JI : in the units of this section, they are 



h 

"rnax — 
Gma 



Uin w U/ (72 - I) 2 + 47 2 2 -72-I 



1/2 



(28) 



LUJ_ 



= 2U 



?n 2 10 (7(72-l) 2 + 47 2 2 - 72 - l) (29) 



with 72 = a 2 n 2 o/a 1 ni and 712 = (a 12 /a 1 )y f n 2 o/n w . To 
estimate £; max and G max , we assume that the density pro- 
file of the initial tpi component has the one-dimensional 
Thomas-Fermi profile n; ni = |V>ini| 2 = (Ml — ^ 2 z 2 /2)/Ui 
with m = (3AaiA r /2v / 26ho) 2/3 - Because half of the ipi 
component is suddenly transferred to ip2, we use the den- 
sity rii n i/2 at z = as an approximation of n^o (i = 1, 2) 
in Eqs. (|28|) and (|29|l . For example, the parameters in 
Fig. Bft>) and (c) yield fc max = 0.205, G max J(J ± = 0.0275, 
and those in Fig. Efb) and (c) yield fc max = 0.441, 
G max /uJ± = 0.592. The quantity 2i? 2 /(27r/fc max ) is ap- 
proximately the number of generated domains in the sim- 
ulation. This quantity equals 3.7 for TV = 5 x 10 3 and 17.3 
for N = 5x 10 4 , in reasonable agreement with the numer- 
ical results. The growth time is determined as 27r/G max , 
which gives 229 (in units of luJ 1 ) for N = 5 x 10 3 and 10.6 
for N — 5 x 10 4 . These times approximate the time scale 
for the first rapid growth of the central density shown 
in Fig. 01a) and Fig. Ufa). At later times, the linear 
analysis is not applicable. 



Analogy to the dynamics of condensates with attractive 
interactions 



We point out that the dynamics described here is anal- 
ogous to the collapse dynamics and soliton-train forma- 
tion in a BEC with attractive interactions 0, fl^L . 
This analogy is as follows. The total density ut = 
l^il 2 + I ^2 1 2 hardly changes during the time evolution 
as seen in Figs. |3fc) andUJc). Thus, we can rewrite the 
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dynamical equations ll'.'il) as 

.dipi ( 1 d 2 A 2 , 1 



Of 



2dz i + Y z2 + 2 {Ul + Ul2)nT 



+l(Ui - U 12 )\^\ 2 - 1(17! - U 12 )\^ 2 \ 2 jr., C-!()a.) 



.dip 2 



1 d 2 X 2 1 



-\{U 2 - U l2 )\^\ 2 + i(C/ 2 - U 12 )\i> 2 \ 2 ),-,. (:50b) 



In this formulation, the term A 2 z 2 /2 + (Ui + Ui 2 )riT /2 = 
Vf (i — 1,2) functions as the nearly static confining 
potential. Then, it is easy to understand how the non- 
linear terms in Eq. I|30() work. If U\ < U\ 2 in Eq. 
(Pnt)j the intracomponent coupling becomes attractive, 
whereas the intercomponent coupling becomes repulsive. 
This favors a spatially localized structure of the i^x com- 
ponent and phase separation between the two compo- 
nents. The same argument applies to the ip 2 component. 
Even when U 2 > Ui 2 , the ip 2 component forms a domain 
structure if Ui < XJ\ 2 is satisfied because the modulation 
develops out-of-phase. 

This interpretation based on the condensates with at- 
tractive interactions can be extended to the effective one- 
component description of the domain formation [2^, • 
Particularly, in the case of Stenger et a/.'s experiments 
[l8j. the two-component condensates of 23 Na atoms are 
characterized by U\ ~ U\ 2 = U, in which Eq. (|3L)|) can 
be reduced to 

Then, Eq. I|31b ) is a linear Schrodingcr equation and an 
effective attractive interactio n ap pears for the i\) 2 com- 
ponent because U 2 — U < [23, |4j] • This means that 
the sudden population transfer from |1) to |2) is formally 
equivalent to the sudden change of the atomic interaction 
of 1 2) from positive to negative. Therefore, the generated 
domains in Fig. [21 and 0] may have a solitary wave struc- 
ture such as a bright soliton train in a single component 
condensate [III Il2l fl3 | . 

Some dynamics are also similar to those in the nu- 
merical simulation of bright soliton formation, in which 
bright solitons are first generated at the edge of an ini- 
tial condensate [Til IT3 . fl3| . Since there is no noise in our 
simulation, the MI is triggered by self-interference fringes 
of the wave function. It was shown that the wavelength 
(amplitude) of self-interference fringes in the initial wave 
function is longer (larger) at the edge of the condensate 
than that at the central part 0, ^| . These fringes can 
be the seed of the modulation, first reaching the unstable 
wavelength 27r/fc max at the edge. 



C. Numerical results for region (b) 

We turn to the dynamics for the combination (b) in 
Table [I] This combination of coupling constants does 
not appear in other systems described by similar model 
equations. The remarkable feature in this case is the 
existence of bright solitons supported by the intercompo- 
nent attraction even if the intracomponent interaction is 
repulsive. Some features such as stability and collisional 
properties of this new soliton were studied recently [43| . 
These studies also discussed the dynamics of soliton for- 
mation from the initial state that causes phase separa- 
tion. 

As in case (a) , after preparing the initial states ipi and 
ij) 2 that has an equivalent distribution with an inverted 
parabola and the normalization condition J dz\ipi\ 2 = 
1/2 (i = 1,2), we change instantaneously a\ 2 from 
zero to a negative value. In these simulations, small 
wave fragments with a large kinetic energy are generated 
when the wave functions undergo self-focusing collapse. 
These waves spread to the edges of the area of numer- 
ical simulations and the reflected waves from the edges 
make the calculation unreliable. To prevent this reflec- 
tion, we added an absorptive potential with the form 
V\ = Vb(l + tanh[(z — z )/£] tanh[(z + zq)/^}), where z 
represents the position of the numerical edge and we set 
Vq = lOOz and £ = 5; V\ can absorb only the waves that 
reach the edge. 

In this case, the dynamical evolution should be similar 
to what is observed in an attractive single-component 
BEC. This is due to the fact that the intercomponent 
attraction favors the spatial overlap of the wave function 
such that |i/'i| 2 — 1^2 1 2 - Then, the coupled GP equations 
are reduced to 



.dii>i 



1 0* 
Ydz 



2 ^2 



X 2 z 



+ {U l + U 12 )\^ l \ 2 )^ »=1,2. 

(32) 

Hence, we can expect that the MI occurs for Ui + Ui 2 < 
(di +012 < 0). Using the results in Fig. the necessary 
condition for the MI is given by a 12 < a\ 2 = —y/a\a 2 = 
—5.65 nm. 



1. Condensate dynamics in a harmonic potential 

The main feature of the attractive intercomponent in- 
teraction is to make the condensate self-focus on the cen- 
ter of the harmonic trap. FigureEJa) shows the time evo- 
lution of the central density |^i(0,t)| 2 for N = 5 x 10 3 
and several values of a\ 2 . Due to the mutual attraction 
and the presence of the harmonic trap, the overall density 
contracts at the center. Subsequently, the kinetic-energy 
cost of this focusing makes the condensates expand with 
two components repeating this contraction and expan- 
sion. Above ai2 — —5.7 nm, although the central density 
of the condensates undergoes a large amplitude oscilla- 
tion, no spatial fragmentation occurs. Below a\ 2 ~ —5.7 
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FIG. 5: (Color online.) Numerical results of region (b) for 
N = 5x 10 3 , where the condensates are trapped by a harmonic 
potential, (a) Time development of condensate density IV^I 2 
at z = for ai2 = —5.6, —5.8, and —6.0 nm. The insets show 
the details of the evolution during the first and third focusing. 
The development of ijj2 component is similar to what is seen 
in (a) because the modulation develops in-phase. (b) The 
contour plots of the density of both components with respect 
to time ([0, 750], horizontal axis) and z ([-100, 100], vertical 
axis) for ai2 = —6.0 nm. (c) The density profiles of \ipi\ 2 
(solid-curve) and \tp2\ 2 (dashed-curve) for a\2 = —6.0 nm at 
t = 150, t = 375, and t = 600. The unit of time is wj 1 . 

nm the central density collapses into some pulsed wave 
packets, or bright vector solitons, characterized by a sech- 
type form of both components .43] ■ The existence of 
these solitons are ensured by the intercomponent attrac- 
tive interaction, because the intracomponent interaction 
is repulsive. The second self-focusing at t ~ 220 gener- 
ates three solitary waves [see Fig. EJb)]. The soliton at 
the center does not move whereas the other small two 
solitons propagate outward and come back to the cen- 
ter because of the trapping potential. Then, the two 
propagating solitons merge with the central soliton and 
this merging generates a few bright solitons again. This 
nearly recurrent process repeats several times, creating a 
single soliton at the center with the help of the inelastic 
particle loss. 

2. Condensate dynamics in an expulsive potential 

It is not easy to compare the above numerical results 
with the MI analysis in Sec. Ill CI The density inho- 
mogeneity in the numerical simulation has a significant 
effect on the soliton formation, the MI analysis for the 
homogeneous condensate cannot be applicable. Also, the 
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FIG. 6: (Color online.) Condensate dynamics in an expulsive 
potential, (a) The contour plots of the density of both com- 
ponents with respect to time ([0, 600], horizontal axis) and z 
([-200, 200], vertical axis) for N = 5 x 10 4 and a 12 = -6.0 
nm. (c) The density profiles of \ipi\ 2 (solid-curve) and \ip2 1 2 
(dashed-curve) at t = 150, t = 300 and t = 450. The unit of 
time is ujT 1 . 



large particle loss of the first collapse makes the use of 
the particle number N nontrivial to estimate fc max . To 
prevent the focusing collapse, we ran a similar simula- 
tion but with a trapping potential with negative curva- 
ture, also known as an expulsive potential IT^ ] . Fig- 
ure El shows the resulting dynamics for a weak expul- 
sive potential V = — (0.1A) 2 z 2 /2, where we used the pa- 
rameters N — 5 x 10 4 and an — —6.0. In this case, 
a self-focusing collapse does not occur, but the density 
modulation grows spontaneously from the edges of the 
condensate, forming a bright soliton train. This result 
agrees with the single-component result in Ref. 0] . 

The MI condition in this case is given by the same an- 
alytic form in region (a). For example, the fastest growth 
mode is given by Eq. (|28|l and the characteristic length 
scale 27r/A: max = 14.2 for N — 5 x 10 4 and a\2 = —6.0 nm 
is in reasonable agreement with the wave length of the 
growing modulation and the size of the bright solitons 
[see Fig. EIb)] . 

D. Numerical results for region (c) 

We now address the situation in which the intracom- 
ponent interaction in one component is repulsive, while 
that of the other is attractive. A single-component con- 
densate with attractive interactions can generate bright 
solitary waves because of its nonlinear self-focusing ef- 
fect [ll|,[H,[l3- Here we ask the related question: in the 
presence of two components, how are the self-focusing 
collapse and the formation process of bright solitons af- 
fected by the intercomponent interaction? 

The conditions for the numerical simulation are those 
for region (b). To prevent the reflection of the waves at 
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FIG. 7: (Color online.) Numerical results of region (c) for 
N = 5x 10 3 , where the condensates are trapped by a harmonic 
potential, (a) Time development of condensate density IV^I 2 
(upper panel) and \ip2\ 2 (lower panel) at z — for ai2 = 
nm and 0.5 nm. (b) The contour plots of the density of both 
components with respect to time ([0, 600], horizontal axis) 
and z ([-100, 100], vertical axis) for a\2 = 0.5 nm. (c) The 
density profiles of |V>i| 2 (solid-curve), IV^I 2 (dashed-curve) for 
ai2 = 0.5 nm at t = 78, t = 120, and t = 300. The unit of 
time is U3~7 l . 



contraction and expansion as we found previously (Fig. 
EJ. However, the ip2 component forms no spatial pat- 
tern, probably because the particle number is not large 
enough to cause the instability. However, a\% has a larger 
influence on the MI as it increases, eventually leading to 
the spatial pattern shown in Fig. Efb) and (c). This re- 
sult shows the increase of the MI strength through the 
presence of the other component, as found in Sec. Ill CI 

The numerical simulation reveals how the spatial pat- 
tern forms from the increased MI. Initially, the attractive 
ip2 component focuses at the center. Because of the inter- 
component repulsion, the density modulation grows out- 
of-phase between the two components and the focused ip2 
creates a density dip in tpi at the center. This strong den- 
sity disturbance generates counter-propagating density- 
kinks, also known as dark solitons, in the ip\ component. 
A similar formation mechanism for the single-component 
system was found in Ref. |44| , where the disturbance was 
given by an external potential. The rigidity of the dark 
solitons is ensured by the fact that they are clearly visible 
and propagate stably in the subsequent time evolution, 
as seen in Fig. E[b). The bright solitons of the attrac- 
tive ip2 component also co-propagate, being embedded by 
these dark solitons. This composite soliton is referred to 
as "dark-bright soliton" or "gray-bright soliton" , which is 
characteristic of the system having a vector order param- 
eter [32, EH • Though the composite solitons propagate 
outward, the bright solitons slip out of the density dips 
of the dark solitons, coming back first to the center. This 
causes the collision of the multiple bright solitons, which 
generates again the new composite solitons. A further 
increase in N or 0,12 increases the number of generating 
solitons and their collision gives rise to more fine-density 
ripples. 



the numerical edge, we used the absorbing potential Vj. 
The initial states ipi and tp2 were distributed equally with 
Ni = N 2 = N/2 (i.e., J dz\^\ 2 = 1/2 (i = 1,2)). We 
changed the value of a 2 to a negative one a 2 = —0.2 nm 
at t = . Then, the ip 2 component generates bright soli- 
tons via MI induced by the intracomponent attraction. 
If <Zi2 = 0, the time evolution of ip 2 is the same with 
the single component problem. However, the presence of 
the second component and the resulting intercomponent 
interaction changes the dynamics significantly. In this 
subsection, we consider the case in which the intercom- 
ponent interaction is repulsive. 



1. Condensate dynamics in a harmonic potential 

We first ran the simulation in the presence of a har- 
monic potential for the particle number N = 5 x 10 3 . 
The time evolution of the central density is shown in 
Fig. Ufa). For a\ 2 — 0, while the repulsive ipi com- 
ponent makes a breathing oscillation caused by a sudden 
population change at t — 0, the ip2 component undergoes 



2. Condensate dynamics in different trapping potentials 

Because the focus here is on nonlinear dynamics in- 
duced by the MI from initially miscible condensates, it 
is desirable that the two components are overlapped as 
much as possible while the MI occurs. To prevent the 
focusing collapse of the attractive component, we can 
also consider an expulsive potential. However, the use of 
the same expulsive potential for both components has a 
negative influence on the repulsive component because 
this component quickly expands and disappears. To 
avoid this problem, we use different trapping potentials 
\\ = X 2 z 2 /2 and V2 = \ 2 z 2 /2 for the two components. 
This situation can be realized experimentally using the 
difference in the magnetic <?-factor or the index of hy- 
perfine sublevels for atoms in each component [3(j • Also, 
tuning the wavelength of an optical laser beam can create 
an optical pot ential that depends on the atomic hyperfine 
spin state | 4q . 

To better understand the role of the ip 2 trapping fre- 
quency on the dynamics, we ran a simulation with the 
parameters N — 5 x 10 4 and Q42 = 0.5 nm and the same 
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FIG. 8: (Color online.) Contour plots of the density of the 
attractive component with respect to time ([0, 500], hori- 
zontal axis) and z ([-200, 200], vertical axis) for N — 5 x 10 4 , 
ai2 = 0.5 nm, and (a) A2 = 0, (b) A2 = 0.3A, and (c) 
A 2 = 0.6A (A = 0.02). (d) The density profiles of |V>i| 2 
(solid-curve), \4>2\ 2 (dashed-curve) for the plot (b) at t — 150, 
t = 300, and t = 450. The unit of time is loJ 1 . 
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trapping frequency Ai = A for tpi, but varied the trap- 
ping frequency for ip2- The results are shown in Fig. [S] 
In (a), the time development of only the attractive ip 2 
component is shown for the case without a trapping po- 
tential A2 = 0. Because of the repulsive intercomponent 
interaction the ip2 component is pushed aside by the ipi 
component. In addition, although the modulation grows 
from the edge, the bright solitons cannot move inside so 
they instead go outside. As a result, to cancel the effect 
of the intercomponent repulsion and to make the two 
components overlap, we had to use a trap with a weakly 
positive curvature for ip 2 - We found that for A 2 = 0.3A, 
the trapping potential A|z 2 /2 and the intercomponent 
repulsion E/^lV'il 2 — C^i2^mi/2 are balanced, creating a 
flat effective potential. Then, focusing of V>2 at the cen- 
ter does not occur and we can thus study the pattern 
forming dynamics. Further increase in A2 focuses the tp2 
component into the center as shown in Fig. [HJc) . 

From Fig. Efd), we find that the generated solitons also 
have the gray-bright character, where the bright solitons 
of ip2 combine with the density dips of ipi . The formation 
dynamics is similar to those found in the previous sec- 
tions, in which the out-of-phase modulation grows from 
the condensate edges and evolves into the solitary waves. 
We find that the Mi-induced dynamics is very sensitive 
to the change of 012- For ai2 7^ both the characteristic 
time scale for the MI to start and the wave length of the 
initially developed modulation from the condensate edge 
[the left panel of Fig[HJd)] decrease from those found for 
the simulation of 012 = (not shown). As a result, the 
size of the bright solitons become smaller with increasing 
a\2- This is the multiplication effect of the MI caused by 
the intercomponent interaction. 



FIG. 9: (Color online.) Numerical results of region (d) for 
N = 5x 10 3 , where the condensates are trapped by a harmonic 
potential, (a) Time development of condensate density |V'i| 2 
(upper panel) and \ip2^ (lower panel) at z = for N = 5 x 10 3 
and 012 = nm and —0.5 nm. (b) Contour plots of the density 
of both components with respect to time ([0, 600], horizontal 
axis) and z ([-100, 100], vertical axis) for 012 = —0.5 nm. (c) 
The density profiles of \^>i\ 2 (solid-curve) and 1^2 1 2 (dashed- 
curve) for 012 = —0.5 nm at t = 80, t = 120, and t = 300. 
The unit of time is loT 1 . 



E. Numerical results for region (d) 

Finally, we discuss the dynamics when both 0,2 and 
012 are negative. The numerical procedure is the same 
as that in the last section except this case has different 
scattering lengths. At t = we give a 2 — —0.2 nm and 
some negative value of ayi- 

1. Condensate dynamics in a harmonic potential 

According to the MI analysis, the modulation should 
develop in-phase in this region. Thus, in the harmonic 
potential, the density focusing due to the attractive com- 
ponent creates a local density hump of the repulsive 
component. Both sides of this density hump evolve to 
a counter-propagating dark soliton pair in the repulsive 
component, as seen in Fig. Elb). Before these dark soli- 
tons come back to the center, the second focusing collapse 
occurs and generates a new dark soliton pair. On the 
other hand, the focusing collapse also generates counter- 
propagating bright solitons in the attractive component. 
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nent. Then, a train of composite solitons forms through 
the MI. Each soliton is largely the bright component of ip2 
with a small fraction of the component being trapped 
by the bright soliton. The trapped if)\ component creates 
the density peaks upon the bright soliton despite the re- 
pulsive interaction and the peaks persist during the free 
expansion. This soliton has been called a bright-antidark 
solitons 0- With increasing |ai 2 1 , the expansion be- 
comes slower and a larger number of sharp composite 
solitons are formed. A similar trapping mechanism of 
the solitons is seen in a mixture of bosons and fermions 
|47|: however, the density expansion of the fermions in 
that case is caused by the Pauli exclusion principle. 



FIG. 10: (Color online.) Simulation for N = 5 x 10 4 and 
di2 = —1.0 nm, in the case without a trapping potential, (a) 
The contour plots of the density of both components with 
respect to time ([0, 500], horizontal axis) and z ([-200, 200], 
vertical axis), (c) The density profiles of l^i] 2 (solid-curve) 
and Iffel 2 (dashed-curve) at t = 140, t = 250, and t = 375. 
The unit of time is luT 1 . 



This combined dynamical process and the resulting mul- 
tiple collisions of the solitons at the center make the dy- 
namics extremely complex. As the particle number is in- 
creased, although the above dynamical feature seems to 
be similar, we cannot obtain a clear physical picture and 
the reliable numerical accuracy because of the rapid den- 
sity fluctuations generated through the above processes. 

As seen in Fig. EJa), the density focusing at the cen- 
ter occurs faster than that for = 0, while in region 
(c) it occurs slower than that for a\% = [see FigEfa)]. 
This difference is caused by the effect of the intercom- 
ponent interaction; the mutual attraction quickens the 
focusing process of the component, whereas the mu- 
tual repulsion delays the focusing. Another feature is 
that the dynamic behavior of the dark solitons in tpi are 
independent of the behavior of the bright solitons in ip2- 
These solitons do not form composite vector solitons, be- 
cause the mutual attraction acts against the coupling of 
the density dips and the density peaks. 

2. Condensate dynamics after turning off the trapping 
potentials 

To see the Mi-induced dynamics more clearly, we con- 
sider a simple situation in which the two components 
freely expand by turning off the trapping potential. We 
first prepared the same initial condition as in the pre- 
vious section (jlll E 1(1 and then turned off the trapping 
potential at t = 0. Figure ITUI represents the results for 
N = 5 x 10 4 and a 12 = -1.0 nm. Because of the mutual 
attraction, the repulsive component expands more slowly 
than that with the simulation with a\2 — 0. After that, 
the MI starts from the edge of the attractive component 
and the modulation becomes in-phase for each compo- 



IV. CONCLUSION 

We analyzed the modulation instability (MI) and the 
nonlinear dynamics of multiple solitary-wave formation 
in trapped two-component BECs. The MI of this system 
was classified according to the signs and magnitudes of 
the s-wave scattering lengths. Then, we used the one- 
dimensional coupled GP equations with the particle loss 
term to numerically simulate the nonlinear dynamical 
stage after the MI occurs. For each combination of the 
scattering lengths, an unstable modulation grew up to 
a train of vector solitons unique to the two component 
system. As a result, we obtained the following picture: 

(a) When all coupling constants were repulsive, the 
strong intercomponent interaction caused the phase sep- 
aration of the two components. The MI first grew near 
the edge of the condensate, giving rise to solitary waves 
with alternating condensate domains. These phenomena 
reproduced the experimental observation by Miesner et 
al. ■ Because the density modulation developed out- 
of-phase, the total density hardly changed during the do- 
main formation, and this allowed us to reduce the system 
to a single-component condensate with attractive inter- 
actions. 

(b) When the coupling resulted in strongly attractive 
intercomponent forces, the two components underwent a 
focusing collapse despite of their repulsive intracompo- 
nent interactions. If the condensates were moved to an 
expulsive potential, the instability of the in-phase mod- 
ulation generated a vector bright soliton train. In this 
case, the fact that the two components always overlapped 
reduced the system to that of a single-component conden- 
sate. Thus, the dynamical feature was similar to that in 
a single-component case [l3|. 

(c) When one of the components had an intracompo- 
nent attractive interaction, the presence of the other com- 
ponent increased the growth rate of the MI over that of 
the single-component case. Moreover, the unstable dy- 
namics was sensitive to the shape of the trapping poten- 
tial. For a harmonic potential, the density of the attrac- 
tive component focused in one spatial region, and this 
focusing strongly perturbed the repulsive component in 
this region, which subsequently produced a train of dark 
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solitons. Then, the attractive component coupled with 
the dark solitons such that the total condensates formed 
dark-bright solitons. When the trapping potential was 
different from each other, the MI occurred from the den- 
sity edge of an initially miscible condensate, leading to 
the formation of a train of dark-bright solitons. In this 
case, the formation dynamics was greatly influenced by 
the intercomponent repulsion. In particular, an increase 
in <2i2 increased the strength of the MI and increased the 
number of the solitons. 

(d) When the intercomponent interaction was attrac- 
tive, the dynamics became more complex. In a harmonic 
potential, the focusing collapse generated dark solitons in 
the repulsive component and bright solitons in the attrac- 
tive one. However, these solitons could not be coupled 
through the intercomponent attraction. When the trap- 
ping potential was turned off, the in-phase modulation 
developed composite solitons in which some fractions of 
the repulsive component were trapped by the bright soli- 
tons in the attractive component. 

Finally we discuss the feasibility of observing the above 
results experimentally. The situation is that in which 
two components initially have the same density profile 
with an inverted parabola and they also have the same 



position. This can be achieved experimentally by using 
atoms of the same species but different hyperfine levels. 
Then, with an rf-pulse, one can instantaneously transfer 
half of the condensed atoms in one hyperfine level to the 
other level |la. 29], In addition, the Feshbach resonance 
during atomic collisions depends on both the hyperfine 
level and the magnetic field. Thus, a suitable choice of 
the atomic hyperfine levels and control of a magnetic field 
can realize the parameter regime in Table [I] 

Another way to control the MI condition for two- 
component BECs is to use a periodic potential |48| , that 
can change the effective atomic mass. In particular, if 
the atomic interaction is repulsive, the negative effective 
mass corresponds to the anomalous diffraction regime 
and, as a result, the system is modulationally unstable. 
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